{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5,
     "status": "ok",
     "timestamp": 1649954662434,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "5OzU9RtB9fWZ",
    "outputId": "146722e2-c641-46dc-a690-01e8eed9160c"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "# 定义状态转移概率矩阵P，6x6表示有 6 个状态\n",
    "\n",
    "P = [\n",
    "    [0.9, 0.1, 0.0, 0.0, 0.0, 0.0],\n",
    "    [0.5, 0.0, 0.5, 0.0, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0, 0.6, 0.0, 0.4],\n",
    "    [0.0, 0.0, 0.0, 0.0, 0.3, 0.7],\n",
    "    [0.0, 0.2, 0.3, 0.5, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "]\n",
    "P = np.array(P)\n",
    "\n",
    "rewards = [-1, -2, -2, 10, 1, 0]  # 定义奖励函数\n",
    "gamma = 0.5  # 定义折扣因子\n",
    "\n",
    "\n",
    "# 给定一条序列,计算从某个索引（起始状态）开始到序列最后（终止状态）得到的回报，chain[start_index--len(chain)]状态开始计算到序列结束的总折扣回报\n",
    "def compute_return(start_index, chain, gamma):\n",
    "    G = 0\n",
    "    for i in reversed(range(start_index, len(chain))):  #reversed()函数用于反转序列，range(start_index, len(chain))生成一个从start_index到len(chain)-1的整数序列\n",
    "        G = gamma * G + rewards[chain[i] - 1]\n",
    "        # print(f\"计算状态 {chain[i]} 的回报 {G:.5G}\")\n",
    "    return G\n",
    "# 一个状态序列,s1-s2-s3-s6，计算过程类似反向传播\n",
    "chain = [1, 2, 3, 6]\n",
    "start_index = 0\n",
    "G = compute_return(start_index, chain, gamma)   \n",
    "print(\"根据本序列计算得到回报为：%s。\" % G)\n",
    "\n",
    "# 根据本序列计算得到回报为：-2.5。"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAEvCAYAAABG0bjWAAAgAElEQVR4nO3dbWwb150u8IdtNtv0JRalfLmGkXU4DNAiyZWbSE5hkEaibjWsuxvYiFNRcdBobwI7ZLPBNjXEWla6m4XsdFgnfUFCyUKMqkBiMogKF2nlisqtnTVZo5HlvdJNArcoh/Uahu8X0ZaSdtNm2z33g3omM+QMRVKUKI2eH0DYmjeeGQ7nf97pEUIIEBERket8pNEJICIiouXBIE9ERORSDPJEREQuxSBPRETkUgzyRERELrViQT4ajcLj8Tiu93g8iEajS3oPj8eDVCq1pGMQERG5xaJBPp/Plw2efr8ffr9/0WMMDg5C0zTb9fF4HACQSCQWS44jeYxwOFzzMSp9n1AotKRjhEKhmjI00WjU8b09Hk9F6cpms/B4PMjn8yXr4vH4op9lvaVSKcfMX7l1RES0uEWD/JEjR6Aoim3wjMfj0HUduq4jm82WPQYAxGIxeDweyysejyMWiwFAyToZuGRgsntJw8PDjscol7ZapNPpigO9XVDt6OjA4OBgSaCVQc0uAAPAxMQEOjo6SpbLDEM6nV70XM+ePQtFUeDz+SpKv51QKOT4edi9ijOIoVDIyJSZyQylPIdLly5BVdWa00lEtO6JMnRdFwBEJpOxXQ9AJJNJoaqqcDpUJpMxttM0TWiaZlkfiUSMfVVVtX0veQxd14WqqkLTNGOZEEIkk0ljffGxFUUpd4pl01zrq/gc5PmbKYoiIpGIZZmmaYtex+JzlOeeyWSEpmmLnq+iKELTNOO6V/oqTv9i5D2hqmrJOl3XhaIoIplMWj47+T7yHFVVrfp9iYjoQ2WDfCQSsX1IC7HwADavUxTFNsAAMIJZcaah+G+nQLZYkFcUpSSgyPVOGZRynNKxGKdMkaqqJQHdLmNSfE3NIpFIyTHMGahKjiHT53RelWQSFiM/i0XyjwZ5HZzSavcqvg5ERGTPsbo+m81icHDQaCePRqNGe208Hkc6ncb4+LixfS6XAwBLFXo8HoeqqsYxfD4fkskkBgYGAACjo6NQFAWBQAAAEAgEoGkaRkdHnZJVQlb7CiGQTqeNv4PBIDRNM47dSLJ63iwcDkNVVUxOThrL0um0bXU8AAwODmLz5s3G39ls1jhHc1PK+Pg4crmcbXPC6OgoVFWtuare3PQQjUZLquH9fj90XUcmk4H4y2zJfr/fsp9skpCv7u5uACip3pdpFQsZUQghoCgKNE1bUt8NIqJ1xSn6AzCq1mVpy1y96lRCltXAmUymbGms0pe51G5Xko9EIkZJ1lzNXtws0EiLNXsI4dzkYF4nz8l8jZ3IEnXxMnMpv9Iqe5km8+cg15lL1bIpQAJgVMvbnZv587JLf3GNhLwHiYioMrYlednxSXaU6+7uRiQSAQB0d3cjmUw6lpATiQSSySSCwSAmJyctJbHil6ZpUBSl7Da9vb1lMymJRMIoyZ49e9ZYHovFSjqwxeNxNDc3l+0kdvToUfj9/qo6ltl1FjTz+XxQVdWowbAzMjLi2CGuv7/f+H80GsXExASEEGVrKXK5HDRNs5SgdV0v2a64tGx+ZTIZx+MLIZBMJjE4OFjSI192oFMUBblcDps2bSrZP5VKIRgMQlEUALB0uJMdOmXtkDwmANtjERGRg8VyAeY2XkVRRCaTMTpVOb0qbcuutA24XEc483rzsWQnNrlNMpkUQ0NDQoiFEuH09LQQQojp6WkBQMzOzlaUZnkdqq0pKNfOb+50ViyZTBodCIvf0+76OV1T2WeiXiV583mpqmp0ppO1LOY0FO9nrokw1wzJa2BX+jf3wSAiosqUHUJX3Paey+UQCAQwPj5etvRdXBqVE+EUv2KxGHRddywRFw8H03UdqqpC0zSjlOnxeBAMBiGEsJT8ent7IYRAJBIx+gns27cPMzMzAIDW1lYAC+3giqKgpaXF9hrIYX5LFQgEoKqqMZzQTC6zG6bY3d2N/fv3L+m95fj3zs7OknW1luQleT+YP/NAIGD5LMxCoZBtTUQgEEAmk0F3dzdyuZxxLa5cuQLgw6F/RERUOccgL8evV/Kgr4RdMHGqrrerVnYi93GSSCQghDCCxoULFyxjr0+dOoW2tjbbfWUV8bZt2ypOTzmJRAKDg4OWzIvs4Gg3UVA2m0UkEinbUa6SyWsuXry4bJ3VZNOGbKI4depU2fkJZMdAYCHgj4yMGJ9fIBCwfJaqqhpNMKdOnbLNpBARkTPHIC8nlwkGg5Zez8W9o+1eS52etlrl0lQcBM+cOYO77rrL+DudTuPOO++0Pa7s+b5x48aSdXYT+yw265zP54Omaejp6TGW9fT0QFEU274HgUCgLsE5kUg4ZhTS6bTjtQsGg47HlMEdWMhoydqejo4OIwMn7x0zc61OOp0ueX/z9evo6DDuw3Q6jQcffHBJ14GIaL1xDPK5XK6khG2uTnaq4q10hrJsNmtU1xdXh8sqWrvgWo5dTUGxiYkJ3HHHHQBgVN07ldTPnDkDALYBUtM02/M3Dyu009vbi87OToRCIYRCIei6jomJiarOU7p48eKSZq4Daquuj0aj0HW9pInETDaXALBkFhKJhFFTI4fayZeiKJYhhLt374au60amcTUMhyQiWksa8it0qVQKAwMDRnU9sNA8IMdTB4PBJY3ndpLP56HrOj7zmc8AWKi6B2D8XUwG33q0yZtt377dKMUmk8maz3NiYsIydn6lyCaQYrqu4+abb7Yss8ssyD4I5kxcNpuFruvYvXu3sUyOSij3uwdERORs0SBvnqe81mr4RCJhlHCj0SjOnDljKfHKqurjx4+XLRErioJ0Oo1YLFa2KtnJsWPHAHzY6e7SpUsAgPn5eRw8eNCyrRxylkwmEYvF6vbrdqFQyBiSGIlE0N3d7ThXfTkyfeagWItaq+uBD4fKyZfTbxwUe/DBB6EoChRFMTJQPT09ZfsfXLx4sfqTIyJa5xzHycsHd39/vxF4l9I+LAPC9u3bbY/T29uL7du3l/2BFllFXGnP7+L3HxwctDQn7N69G16vF52dnfjiF79o2b6/v9+YTS6ZTKK7u3tJvz4Xj8eNdmh5LWXVtaIoVR+7v7+/LrUdS+ld7/P5LNs7Vd0Xk73vdV03+jboul5yX8jRHZqmlXRYJCKiClQ75k6OXy73Kp5b3G72NSHKj+muJmnl0lTLXOxOP2xTzY+6FI8JNy+zY54LwG4mO/M4eTlPgd3xys09IH+LoJLPsPglfzyo2v3srknx9ZSft/m+kfMcyGth/iEjIiKqTM1t8qKCjney9DowMFB2mFsxWcpbrLf6YmmqpR1XlvjtOsPJtuhKXkeOHLH0Phc28weYyeFjQgj09PQ4nns+n0cul1tSW344HK74PORrfHzcksZqXsU1AvK+GBwcNGa2E0JgYmIC8Xgc0WjUGL4pO9slEgmoqlrSW5+IiJx5RDXRl2iJQqEQ0uk0dF23zaTI+RmcbkuZ8VlsFAMRETHIExERuVZDhtARERHR8mOQJyIicikGeSIiIpdikCciInIpBnkiIiKXYpAnIiJyKQZ5IiIil2KQJyIicikGeSIiIpdikCciInIpBnkiIiKXYpAnIiJyKQZ5IiIil2KQJyIicikGeSIiIpdikCciInIpBnkiIiKXYpAnIiJyKQZ5IiIil2KQJyIicikGeSJalQqFAqLRKJqbm+HxeHDw4MFGJ4lozWGQJ6JV6dixYwCAq1evYmxsDIcPH0Y2m21wqojWFo8QQjQ6EW5VKBSwZ88epNPpknWqqmJ8fLwBqSJamzweD6anp9Ha2tropBCtGSzJL6PXX38diUQCqqpC13XMzs4CAIQQtgG+UCggHo9XVC0ZDocRj8dRKBTqnu6V5PF4jFc8Hq/5OPl8HtFoFKlUquZjxONxRKNR5PP5mo/RSPW4Bqv1vjp58iQikciSA/x6vE/i8fiSvltrUTabrduzZc0TVDUAxiuTyQghhJidnbUsn52dNbZXFEUIIcTY2JhQVdX2mLOzs6KtrU0kk8mK0zE0NCTa2tos77UaqKpquRbFL03TjG2L/y53jLa2NjE9PV2y7fT0tFAUxXZdtTKZTE3HymQyjp9trRp1DWq9r1RVNb4P9TQ9Pe14bXmfLE7TNNvv2Hrh9IxZL1iSr8HY2BgAQNM0BAIBAEBLSws0TQMAZDIZtLS0AFjIUba1tQEAfvGLX6Cjo8P2mHv27MEDDzyAcDhccTr27duHBx54AF/96ldrPpfl8PLLLwNYqLFQVRWZTAZCCAghEIlEqjpGMpmEEAKzs7NQFAX33nuvZbtCoYB7770XIyMjdanGDQQCGBkZwf3339/w0myjrsFquq9mZmYQi8WMa1GM9wlReQzyNdixYwcA4Pbbb7csv3jxoiXwA8Dx48eNB8758+exbdu2kqqjbDaLyclJ9Pb2Vp2WRx55BK+88sqq6pAkMzh2EolERed54cIFAMDWrVuNY+7cuRPXrl2zbHfs2DF0dnZarvlSBQIBeL1eo+NXozTyGqyG+6pQKODRRx/Fyy+/jJaWFqRSqZL08D4hKo9Bfgnefvtt4//ZbBbnzp2zBLB8Po/BwUF87nOfAwDcdddduO+++3DzzTdbjnP8+HHHEvzBgweNIUQejwfRaNSyvqWlBV1dXTh+/Hi9TguhUMjSnmX3CoVCVR+3vb294m3Pnj0LRVHg8/mMZW+99Ra8Xq9lu+HhYezcubNk/0KhgHA4bKS3ubkZR48erfj9H330UQwPD1e8fTVSqVTJZ5nP50uGiS31GgCL3z9OluO+MotGo/B4PJa2cXldTp48CQB46qmnMDU1hZtuugkejwfd3d0lx3HzfbKcZmZm0N7ebpy33++vKUNXj2dFKpWC3+83tm9vby9bO1IoFBZ9X/oQg3wdFAoF9PT04MUXX7Qs9/l8EEIY1YOHDh3C1atXSwL6xMQENm/eXHLcbDaLw4cPY2pqCkIIDA0N2b7/nXfeiXPnztXpbIDx8XGjet3pVc3IgGAwCI/Hg6mpqYr3efXVV9HZ2Qlg4foePXoUhw8fRiKRMLbJ5/PQdR2bNm0q2f/YsWNGZ0chBDo7O7Fhw4aK3/+2226Drut1r4rN5/OYn59HMpnE4OCgsfzIkSNQVRWHDh0yli31GlR6/zip930lpVIp7Nu3D5FIBCMjI8by/v5+aJpm1JQlEomS+664JO7W+0Rqbm62ZPyOHj2K5uZmS6e/4m0qEYvF0N7ebjRxAMDGjRurTl89nhXd3d0YGBiAEALT09O4du1a2drA5557DolEAn19fRgaGoIQAmNjY1AUxXhPMlneJn/3amtrMzpz9PX1LaljB0wd+Mymp6cFgEU742maJlbbRynTY+6M1dbWZrtd8bUr7sQIwLZTVyaTcTzvoaEhoSiK43VVFEUAEJFIxLaDmTy2U0cyec0Xezkxp1124jKnox7XoNz9Mzs7KzRNE11dXY5pXOy+quT8y30vNE0zOqEtlhY7y32faJpm3CeqqjbkPvF6vaKvr8+SXq/XK3Rdd9zGLg3Fn0NXV5dQVdVyHKmvr08AEIqiiLGxMcfj1otMf7UdPc0dLCORiIhEIrbbLXYfut3qigxriKqqQtM0kclkbINXNco9JMbGxoTX63XsMSxE/YP8Yr3j5UOvHLsg77Rd8RdwbGxMADDOt6ury/Yal3t4C/Hhdenq6rI8QLxerxgbGxOzs7NCURQxNDTkeOxKe4tX22taBighhO1nW69rYHf/TE9PGw/Fcmmu9r6qtne9HG2i63pNvfmX8z6R3+vZ2Vmh67rwer22gWK575N6sAvys7Ozoqurq+T7NzQ0ZGRoZIainHo8K6anp0VbW5vxvayE+fsjxMIIJqfC0HoP8qyuX4L5+XnbavpavPvuu7bLd+zYgatXr0JRFDz66KNLfp9K1Lu63iyVSi06ZvWnP/0pFEUxmjkef/xxTE1NYWZmxnZ7p6rS3t5ezM7OYmpqCs899xwAGO2OO3bsQEtLC/bu3YsTJ07UdC5LIasj4/E4HnjggZIe3/W6Bnb3T2trKxKJBB588MF6nU5NbrzxRgAL7fPf+c53ylbR2lnO+yQQCODcuXNoaWmBz+fD1q1bq6rGX+1kR8bp6WnEYjGjH8SJEyewa9cutLS0YN++fbh27VrZtvp6PCtaW1tx7tw5PPPMM/jSl75UUdPH66+/DlVVAXzYHCM7X5IVg/wSHD58GHv37l3ykBxFUSyd+ICFYGhuu29ubi57DDlMbzVw+pIWCgU8++yzi+4rO+JIsh32l7/8pWVb2YYoe1hL0WjU0nnK3AnrnXfesTwMikdIFPvMZz5Tdn2t5DU6depUyWiDelyDau8fJ8t5X7377rtIp9Po6Oioutf7ct8nZrL9e/fu3Y7pWa77ZDn4/X4jI/SpT33Ksm5ychK33Xab8fdyfv75fB5+v9/4LthlopwmshkZGTGGI1+5csXYn79vUIpBvkYdHR1oa2uradhbsc7OTly8eLFk+dTUlNFbNJ/PO9YYnDp1qqqe68ttz549ABa+oOl02uh4d9NNNy3a+W7Pnj24du0a0um05cvd1dWFAwcOWHpj+3w+KIqCy5cvW47h9Xpx4MAB4z3b29vx5JNPAliofTGTpclistd2taXLSl24cAFer9fSQUyqxzUAKr9/nCz3ffX222/X/B1a7vtEkj+Sk0gkbO+F5b5PloPf78eWLVvg8XjQ1tZm6exYPPRwOc/rvffeg9frNUZP9Pf3Y2xsrOQ9i4P/yZMnkU6njZL8xo0b4fV6cffdd+PLX/7ysqV3zVr5FgIqlslkhNfrrWnmOtk2tRwzja0ELLG9zNx5qxLyWi+2v7ljZaXHrSYdyWSy6o5mTqq9BlK5NNdyX1XbJh+JRGz7QyyHWq7R7OysUFW17Kx2y32f1EM1M96pqmr5TBr5bJmdnS3pZFiLpT5j1joG+VVCduSrVl9fX92CRSMs9QsoHwSVPojk9slk0nEq4bGxsZLe7vWmKErdHjzVXgNJ9h63O8/lvq90XV/RAFLLNTIHeF3XbTuILvd9Ug/VBHlN04wOh8lk0piSuxH6+vrqMgUxgzytCrXMXZ9MJlfl3PXVqMcXsNo5yZPJpPB6vQJAydCjes5v7kQObatngKs23Sgz1G0l7quhoaEVH/ZZzTVKJpNlr9FK3Cf1Uk2Ql7UX+MsQurVwfothkKdVQ45dLjfmVYpEIkLTtDUd4IUQjg/Raum6LiKRSFWZpGJDQ0MiEoksuXqwUepxDdxyXznhfbI+yKGN9Xi2rHX8PXkiIiKXYu96IiIil2KQJyIicikGeSIiIpdikCciInIpBnkiIiKXYpAnIiJyKQZ5ogYZGBhodBKIyOU4Tp6oQW644QZcvXoVN9xwQ6OTQi7x/vvvo7m5Ge+//36jk0KrBEvyRA3S3Nxc0W9nE1WqUCjU/LPC5E4M8kQNcsstt+DXv/51o5NBLvLrX/8at9xyS6OTQavIdY1OANF69cUvfhF/+7d/2+hkrBnbt2/HmTNnGp2MVY99PciMbfJEtOr953/+Jz75yU/i6aefxlNPPdXo5BCtGayuJ6JV7+Mf/zg+/vGPo6mpqdFJIVpTGOSJaE34yEc+gr/6q79qdDKI1hRW1xPRmnDy5Ens2LGj0ckgWlMY5ImIiFyK1fVEREQuxSBPRETkUgzyRERELsUgT0RE5FIM8kRERC7FIE9ERORSDPJEREQuxSBPRETkUgzyRERELsUgT0RE5FIM8kTUcKlUCs3NzfB4PCWvbDbb6OQRrVkM8kTUUIVCAfPz83jttdcQiUQghEBfXx80TYMQAoFAoGSfbDaLcDiMmZmZmt4zn88jFAohlUotNflEqxqDPBEti3g8bpTGQ6GQsfzgwYPG8oMHD6KlpQX79u3D2bNnsX37dgDAxMQEtm3bZnvcVCqFr33ta3jhhRfQ2tpaU9p8Ph9efvllnDlzBtFotKZjEK0F/BU6Ilo27e3t0HUdv/nNb9DS0gJgoeR+6623orOzEy+88IKxPBQKQdM0bNq0CTfddBPsHk3ZbBb33Xef5XhLFQqFsGvXLuzbt68uxyNaTViSJ6Jl88ADD0BRFEtAnp+fBwBLgM/n85icnERraysuXLgAVVUxMzODkydPWo43MDCAZ555pm4BHgCeeOIJHDhwAIVCoW7HJFotGOSJqGqFQgHt7e1obm62tItHo1E0Nzcjn88by6ampiz7RqNRvPTSS5ZAfeTIEXR2dgIANm7ciMnJScRiMdx9993GNvl8Hul0Gl/4whdK0jMzM4P29najGcDv91fcYW/Hjh0AgNdff72i7YnWEgZ5IqrasWPHMD4+jubmZqTTaQALQXhwcBCvvfYafD6f7X5Hjx5FU1OTEVilRCJhdILz+Xy4evUqxsfHLRmByclJY32xWCyG9vZ2CCEwOzsLYCGzUKmtW7firbfeqnh7orWCQZ6Iqtbb24uWlhb4/X5jWTQaxdDQkKU3/O233278P5/P49vf/jZeeOGFmt7z0qVLUFXVdl1TUxPy+Tzy+TxaWlqQy+WMzEA2m0UoFFq0ZH/+/Pma0kW0mjHIE1HN7rrrLgDAyZMn0dTUVNJ57cYbbzT+H41G8f3vf7+u7enSCy+8gKamJiiKgng8biyPx+P42c9+hlwuV/f3JFoLrmt0Aoho7dqwYQPm5+fxxBNP4M0333TczqmavlpOneNaWlqQSqVw4MABbNmyBbfffjt27NiB3t5eACyl0/rFIE9ENZufn8fhw4eRyWTKltAPHDiA3/zmN0t6r5tvvrmkEx8A+P1+/OhHP0Jrays+9alP1Xx8WStB5Casrieimp0/fx59fX22s9IBMJYX96avxdatWwHA0nMfWAjyW7ZsgcfjQVtbGzRNq7rGIJ1O44477lhS+ohWI5bkiahmuVwOmqaV3aZe8235fD6oqorR0VGjGh4AxsfHl3TckydPwuv12g7NI1rrWJInopqkUinoul7z1LK16O/vx7e+9a2qJq4pFAqYnJzEO++8Y7v+n//5n+s+wQ7RasFpbYmoJuFwGHNzc0suSVcrlUrh2WefLRlHbycejyMWi1mWmR95ct76RCJR/4QSrQIM8kS05mSzWTz//PM4cOBATTUJhUIBe/bsQU9PD8Lh8DKkkGh1YJAnIiJyKbbJExERuRSDPBERkUsxyBMREbkUgzwREZFLMcgTERG5FIM8ERGRSzHIExERuRSDPBERkUsxyBMREbkUgzwREZFLMcgTERG5FIM8ERGRSzHIExERuRSDPBERkUsxyBMREbkUgzwREZFLMcgT1ai5uRkej6cuLyKi5cAgT1SjRCJh/F/TNAghqnolk8kGpp6I1gOPEEI0OhFEa1U4HMYrr7wCAMhkMggEAlXtH41GMTg4CH4NiWg5sCRPtAQvvPACFEUBAPT09KBQKFS1//bt25cjWUREABjkiZakpaUFP/rRjwAAuq7jq1/9alX7b9q0aTmSRUQEgEGeaMlaW1sxNDQEAHjllVdw9OjRivcNBAKsqieiZcM2eaI6ke3zXq8Xp0+fRmtra6OTRETrHIM8UZ0UCgXceuutuHbtGtra2jA+Po6WlpZGJ4uI1jFW1xPVSUtLC1577TUAwNTUFJ577rkGp2j1OXjw4JqaF6C9vR2hUKjRyWiId5/4B/zum183/v7dN7+Od5/4hwamyNn8w7vwwel0yfI//eodzD+8C++/+HzVx5T7fnA6jT+Mvmy5FmsJgzw1RD6fRzQaRSqVWnS7UCi06HarRSAQgKZpAIDDhw/j5MmTDU7R8kulUhUHwvPnz6OrqwvAQvNGPB6vekSCWaFQQCgUsp1gyC5Nld530tTUFHbt2oWTJ08iHA4jm83WnNaV8Psj/4r5h3eVvP70q3canbRlIwP49feqJev++NMfwbOhCTc8+viS3uO627cAwKrN5JQliOpAVVUBwPJqa2sT09PTJdtOT08LRVFs19mZnZ0VkUhERCKRqtPVqFtcXg+v1ytmZ2eXdKxMJlNybb1er+jr66tTamtX7ecCQCSTSePvoaEh0dbWVnKNzOeayWSEEAv3gXn57OysSCaTQtd1oaqq0HXd2MZOtfedvO66rgshhNB1XbS1tVnSv9r87ttPi7mv7LRd9l8X3q7qWPP/2CPee+pJ4+/3nnpSzP9jT13SWS9/PDUu5r6y0/Y1/489ZddVYu4rO8UfT40bf//u209b/l4LGOSpLuTDVT4AZ2dnRVdXl/B6vSXbeb1e48FdDVVVxdDQUFX7NCrI67ouvF6vACBUVV3y8fr6+oSiKMbfQ0NDJQFzpWmaVtW5TU9Pl9wP8jhdXV2WZWNjYwKA0DStZFtz4JfktRkbG7NNUy333dDQUEm65HEqzSisNLsg/18X3hZzX9kp3n/1paqOtRaC/NxXdor3nnrSOEez9556Uvzu208LIRbOpTg4l8sELPaSx10LrluxKgNytQsXLgAAtm7dCmChfXrnzp3GbHDSsWPH0NnZWfXMcADwxBNP4KGHHsLu3btXfYc2n8+Hl156CV/60peQTqcRj8fR29tb8/HOnz+Pzs5O4+99+/bhsccew6VLl+qR3KoVCgV861vfwunTpyve58KFCwiHwyXLH3nkEcRiMTz++OPGfbFjxw4AwO23327Z9uLFi9A0zXL/ZLNZtLW1AQB+8YtfoKOjo+Q9arnvZmZmsHPnTsuylpYWRCIRxGIxjI+PV3ysRvrv/3fZdvkfRl/GH38yavx9ffDzVVVr//7Iv+JPb/0f3NDzGN4fGTKW//Xf78bHdu/Bu0/8A8T8HADgo3/jwyf/9VnL/ub15v2Ahfbw3z/Tj+uDn8d/z13Fn976P8Z2G354AsBCH4Hr7vgsPrH/mwCATxwYwPzDu7DhhyeMKny57sbv/wDzD++yVOnf+P0fLHqO7z7xD7h++98a6VqL2Ca/zmWz2Yp+QGWxtsizZ89CURT4fD5j2aI67u4AACAASURBVFtvvQWv12vZbnh4uOTBCSwEjXA4bLxfc3NzyXhz+eB//fXXaz1dR4VCAe3t7WhubsbMzIyxPBqNorm5Gfl8vupj7tixA319fQCAWCxmOW610uk0/u7v/s74W34exUFwpbz++utobm62HSaYSqXg9/uNz7K9vR2FQgGXLl3Cgw8+WLJ9S0sLurq6cPz48ZJ1b7/9tvH/bDaLc+fOlWSWjh8/jnvvvRfAQmZo27ZtiMfjlm2c7jtgoTOg+ceGotEoAODq1av4whe+ULL9l7/8ZaTT6ZruiUb4w4kUPvo3Pkugev/F5/HHn4zihp7HsOGHJ/CJAwP4IPPzmjqofXB6Aht+eAIbfngCH/0bH/74k1HMP7wLH9sVxoYfnsANPY/hz/+Rtxz790f+FR+P7jf2uz74efzxJ6MlfQc+yPwcH92sGNsBH7aLy0yD7Hfw+2f6jb8/yPwcf/6PvKVfgnnb9YQl+XWuXpOxvPrqq0ZJs1AoYHR0FIcPH7b8CEs+n4eu67azvB07dgy6rmN2dhYtLS0Ih8PYsGFDyXZbt27FW2+9ZVsiXIpjx45hfHwcd999N9LpNFpbW5HP5zE4OIhMJmPJvFTj0KFDmJiYwNTUlHHcasnOe3fffTeAhWD3ta99DW1tbUbGZ6WdOXMGfr/fdl13dzeSySTC4TBmZmZw//33o6WlpWxNxp133olXX33VcX2hUEBPT48xu6AkP6Pp6WkAwF133YX77rvP8uNB5e67bDaLw4cPQ9d1+Hw+HD161MiMOXXOk5/hr371q5rvi+VWHMg+WVRq/SDzc1wf/LxRsr3u07fhujs+iw8yP6+6k9rHHvxfxv+vv7cT748M4bo7Pmsc+/p7VXxwegJ/vvRbYztZwpb+KnAvPsj8HP+VPY3rPn2bsfy6Oz5ryZxcH/w8Psj8HH/61Tu47tO3lRxHev/F5/HnS78tqT0w++B02lIDUeyGnsfwkaZmiLlrjtusBSzJ05IVCgVMTU1hcHAQHo8HN910E06cOIFMJmMJxleuXAEA2yrTDRs24Nq1a0a1fyqVcgzk58+ft11uVysBoKJaid7eXrS0tFgCVzQaxdDQUE1NC2adnZ3wer145JFHatr/F7/4BQDgpptugsfjwX333Yf29vYlVxfPzMwgGo06lkiz2SwOHjxouz6fz9tWiwOA1+vFW2+9hUKhgNbWVuRyuYrSMzU1ZflbVsEDwHPPPYe9e/eWZJJ8Ph+EEMbyQ4cO4erVqxXfd5/61KcAAJOTkwAWmkHMGYRyzLUMq40s+W744Qn89d/vxvzDu/CH0ZcBwPj3o8qtln0+unnhNxiq7YlvDsrFxzL777mrtvt/cDptlMKLfaSp2fK3p2mhZtCpCaIa19+rGtfoEwcGACxU+Ru1C/eq8Ny4wTHdawVL8utcNptFMBhcdLtyv7D25ptvAgCmp6fR2tqKcDgMXderCo779u3D/Pw8gsEgurq68MILL1Td7m5XK+HxeKqqqbjrrrsALJSem5qasG/fvqrSUOzkyZM4fPgwMplMzf0IJiYmEIlEKg4+lZiZmcEzzzxTdihZIBDAZz7zGYRCoaom9jl9+jQeffRRDA4O4qWXXqq5tkG+XzabxcTEBM6dO1fTccppbW3F2NgYHnroITz77LN48cUXXTdT4cd278Gf/u+/44Mz/9taZT8yZFuS/e//dxmwCdz19Ltvfh1//o+FzONH/8aHTxwYcAz0lRzDSbmqeVn9X85HmpotNRBrEYP8OleP6vqf/vSnUBTFeDg+/vjjCAaDmJmZsX1gFgoF24DR29uLRx55BHfffTeee+45HDp0aEnpqsWGDRswPz+PJ554wsi81KpQKOChhx5CX19fzbUB+XweU1NT+PrX7SfiKBQK2LNnD9LphYlANE2Dqqq4//77oet6yfaapqG3txf3338/JiYmLOsOHjyIw4cPAwBUVTUC+9e//nU89dRTJZmM+fl52zS1trbi3LlzOHr0KL70pS8ZTTC1mJ+ft62mr4XTfbdjxw6j9P/oo48uS2ai0Tw3boAoCoifODBgWwpfbn8YfRl//o88buh5zKjSr2Ucf7mq+GLzD++yvF+lPE1e/Dnz82qTtqqwup6WpFAoGB2tJBnQfvnLX1q23bhxI4APe+JL0WjU0smuuLNeMVnaXg7z8/M4fPgwRkZGltyDf8+ePVAUBU8++WTNxxgdXej97PRrdaOjo2hqaoIQAmNjYwAWrvubb74JRVGM5aqqQgiB3t5ezMzMwO/3W9qTZ2Zm8Morr2B2dha6rqOpqclYFw6HS0r8Pp+vpNkkn8/D7/cbk9vY9akox1w9Lx0+fNi2mr4aTvddcZNQc7O1angxjer0WAvx7jw8GxY+Uzmxy5/enm5MWv7Sxv2R//HhPV2P6vflIK+V3Wx6awVL8rQke/bswbVr10qGiXV1deHAgQPYsGGD8SD1+XxQFAWXL1u/0F6vFwcOHMBjjz0GAIhEIo6BMZ1Oo6enZ9nO5/z580sqeUtHjx5FOp3G9PR0zZmFbDaLWCwGAAgGg7Y1Lp/73Ofw7W9/G/F4HI888ohRNZ7NZo2OkMXDyi5cuFCSUdq0aRO8Xi+eeuop7N+/vySoK4qCfD5vZAy2b99eUhPw3nvvwev14qabbjL2GRsbq+j8T506hfb2dsuyjo4OFAqFJQ09BJzvO2ChH4Dsu6GqKl588cVFjyf7dHz6059eUrpWiiw5//Xf7wbwYSe7P/5kFB9pucko3f5h9GX8+aLu2JmtXmS7urmTXbkOcE7kMLtKmZsnyg0XNDcDyBED/3X+TeM6VdKpbzVhkKclcer85dTWu3fvXoyMjFhKUIcOHaqoav7kyZPwer22w5rqJZfLGdPS1mpmZgaPPfYYhoaGllQCraQppbW1FRMTE+jr68OpU6eMz+Ps2bPGe09MTOA73/mOsc+lS5dKStktLS0YHx/HsWPH0NbWht/85jeW4NzS0oIrV64YQf4LX/gCotEostmskSGS1fTVKhQKSKfTyGQyluW9vb1LDvCS3X0XDodrGqXxs5/9DKqqrtqe9UBpW3Rx1fwn9n8T77/4vCXwmcecL6eP7d4DMXcNH2R+jg/+UhVePNa+Etd9+jZLu/r8w7scg3e56no53h8Afv9Mf8l21/3PO43hfdd9+jb8+dJv8dGbb6kqrQ3VqFl4aH1ayox3bW1tVc94V41kMrnkGfJmZ2dFW1tbXWa5W0wmkzGmti2eTc48fSv+MgWsvHbJZLJkSlxVVcXs7KyYnZ21nfq1ra3NmN5VqnbGOyd9fX0lM8vV21LuOzM5k+FqnfFuvZJT9r731JOWWfqk4ulphfhwdsBKZgJ876knjW3tjrWaMcjTiqt2DnEhqp8jvRZdXV1LDlp9fX11ma++EmNjY8bUuYqiGFPcJpNJS8Bva2uzXO/p6WnLeeq6LhRFKTsnvt10tEIs/XNJJpO2c9cvh1ruOzOZgVvNc9fTAhnAq52rvhwZ6Ffb1L6LYZCnhtB1XUQikUUfmLOzs0JV1TXxYJXzrY+NjTU6KYtSFKWkZO4kmUyWDeTJZLKmzFEkEhGapq1IgJcqve+KZTIZ0dXVteSaAKKV5hGiDtOdEa1zhUIBt956KyKRSEOG/lUrm83i+eefX/QnV+VPuVYzTp6IVg8GeaI6CIVCyOVyePPNN9dMMMxms/jZz36GRx55xLYT2WLriWj1Y5AnWqKjR4/iscceM2b8q1ZzczNee+21JQ/bIyIqxslwiJZADpfTNK2mAD8zM4Nr19b2D2AQ0erFkjxRjWR7tRxjXsv+X/3qV/HKK68Yv4JGRFRPnAyHVszFixexefPmRiejIpWk9bnnnjN+OU3OmlYrBngiWg4sydOKmJubw2c/+1n89rer/xedKknrzMwMtmzZUrf35NeQiJYD2+RpRXz3u9/FxYsX8S//8i+NTsqiKklra2srxMI8E3V5EREtB5bkadnNzc3hlltuwdzcHJqamlZ1R7OVTms8Hserr76KlpYW5HI5DAwM1DSfOhGRHZbkadl997vfxdzcHICFILqaS/P1TmuhUEA4HIbH4zFeoVAIwMI49OHhYZw7dw7j4+Po6upy/EnZert48eKKvE8jrYdzJFoMgzwtu5mZGdxzzz0AgHvuuQczMzMNTpGzStMqe9bLwB2PxwEs/Pqex+Mx9tuzZw/uvfdeS9W87In/zjvvwOv1Gsc8dOjQioyVn5ubw65du1wdBN944w3s2rVr8Q2JXI5BnpbdiRMncPr0aQDA6dOnceLEiUX2aJxK0zo6OoqmpiYIITA2NmYsP3PmDBRFMcbM53I57Nu3z/YYu3cv/L53KBRCPp+v52mUNTIygunpaTz99NMr9p4r7Xvf+x6mp6fxxhtvNDopRA3FIE9Ug8997nOYmppCPB7H3Xffbfzu+blz5zAwMGBsp+u6pareXF0vx9f7fD60tbWtSKCfm5vD9773PQDAj3/8Y1eW5t944w38+Mc/BgB87Wtfa3BqiBqLQZ6oBq2trZiYmMC///u/Y8+ePQAWqvB1Xbd0nFMUpaQnvXninJaWFiQSCSiKgtHR0WVPtzmwz83NubI0bz6n6elpTE9PNzA1RI3FIE9UpWw2i4MHD8Ln8+HAgQOYnJwEAFy+fBnf+MY3LNv6/X4cPXq05BjxeBwnT560LNu2bdvyJfov/u3f/s0Y33/PPffg4sWLRkdDN5AZmC1btqCpqQlbtmwxSvVE6xGH0NGK8Xg8a2ZMeLm0njx5Eg899BCuXbsGRVHKDnszT10rqaqKu+66C6+88gr8fj8AoKOjw6jyXwm33HILTp8+vWZmIKzWG2+8gaefftroX0G0XjHI04pxS5B3AwZ5ovWB1fVEREQuxSBPtA41NTW5qi2eiOwxyBOR68hpiYnWOwZ5IiIil2KQJyIicikGeaJ1iG3yROsDgzwREZFLMcgTERG5FIM8EbkOe9cTLWCQJ1qH2CZPtD4wyBMREbkUgzwREZFLMcgTERG5FIM80TrENnmi9YFBnoiIyKUY5InIdTiEjmgBgzwREZFLMcgTrUNskydaHxjkiYiIXIpBnoiIyKUY5ImIiFyKQZ5oHdqwYQPb5InWgesanQBaP3p6ehqdBPqLf/qnf2p0EohoBTDI04r5wQ9+0Ogk0F+4fQz53NwcNmzY0OhkEDUcq+uJiIhcikGeiIjIpRjkCQAQCoWQSqUanQxqgFQqhVAo5Lje4/Egn89XfdxQKIRoNLqUpBHREjHIE4CFTnHd3d01Pcxrlc/nSwKI3+9fcmYjm83C7/cvNXkrKpvNwuPxlCyPx+NGAK4m2FaTabt06RI6Ojoc0wUAPp+v7DH8fj88Ho/xnvl8Hul0Gg8++GBFaSCi5cEgTwCAcDgMTdNw5MiRstvJYGT3yufzjuvlwz8ejyMejzseW9d1hMNhx/ePx+OO7+/xeIyg5DbZbBaRSASKolR0jul0uuJjnzp1Ctu2bbNdd/nyZaiqapse83UfGRmBEAJbt24FAKMEHwwGHe8FIlp+DPLrTLkgGYvFMDg46LjeTAhhvHRdL3kf83pFUYzlu3fvRiwWsw30AwMDALBo8FZV1XJ88ysQCNTrUi07+VlUIhAIIJFIQNM04zqV+ywBoLu7u6LPMpfLWYJxNBpFKBSCx+NBd3c30um0Zd9QKIRAIGBcczOfz4dUKoV0Ol3y2SSTSQAom4mrl/n5edePICCqBIP8OtPb2+sYIBd71YPP54MQArFYzLI8m83aBgaZgag0eEejUQSDQei6XlLLsNr09vZCURTHmg2nfcbHx43/231OkUgEqqqWzQwBH2amOjs7LesSiQTGx8eNDFomk7Gsl+9vJ5vNoru7G5qmWZpM5PJMJlPj1SKiWjDIr1PRaHTRdmtZUrQLkOYAai6pV6o40xAMBhGJRErSFI1GoWlayf6yDbj4lUgkjEDiVMuwmuzdu9eS4SkuncdisZKSdLkq73g8jomJCYyPj6Onp6dshzohBFRVxfbt2xEKhUoyG7L2xC6DZU4n8GG1fDAYRDKZxO7du9HZ2WnUwsjla6mmhcgNGOTXqf3790PX9bLtu8PDw1BV1bbT1WLV9YtJpVJG5uDIkSOIRCLYv38/Ojs7jUAv23V7e3tL9s/lcohEItA0re61DStJntvZs2eNv4tL5YqilJTE7aq8/X4/YrEYOjs7ASxUi/f09JRtEkin0yXHikajRsAurhGRr23btkEIYWSozKX9cDgMn8+HRCIBVVWN46xENT0RWTHIr1M+nw+KouD48eO26/P5PHRdr/tUtLLD1pkzZ4zMwf79+5FIJIzAIEvp+Xy+bNWwk8uXL9dUu9AoQgjHjm/ycyjX3CDbzwcGBiCEwMTEhFGCD4fDEELYtsWnUinbTnWJRMKxqr6474PMnDilSZ7fYr3ziWh5MMivY3v37sXg4KDtutHRUQDOnaRqqa7P5/MIBoNGu6+daDSKdDqNZDLpGOBlwMjn84jFYkY67Nq2r1y5YtlntQoGg7bL0+k0VFU1Pg8zmRnq7++3lO5zuRwAWK6HDM7mjowjIyOOmTg50mGx6vXh4WEAC0Mw/X6/Edz7+/uhaRrS6bRrRzwQrQUM8uvY7t27AcC2fXd4eBiRSMRx31qq62Wnu2KTk5NG8Nm+fbtjdbSdZDJptC1Lly5dWlPj5D0ej9Hz3Ew2afT39xvB1CyXyzmOKBgfH7dt5jCXxNPptNEDP51OGxmmaDSKy5cvG2krfsnMQyqVMq7zwMAA/H4/Ojo6jM9Y9jUoHka3Ep/N3Nwce9cTgUF+XfP5fFBVFSMjI5blqVQKuq5XNJFJKBTClStXaqqSlUE/HA4bbbuLBfeLFy8a/5cl1mKnTp1ynNxltfF4PMhkMti0aVPJupGREezdu9cI4naZsVQqVVINX24uAjNzRk1VVaN/QyKRMKr57YbAycyhuSZg06ZNRsZC1thEIpGSY6iqavQZIKLlxyC/znV0dJRMnDIyMgJVVSvqCZ1IJBAMBiseouY0bltWVzutl6W/fD6P7du3G23VcvIVALj55puNmdac2rhXk1Ao5NjjXI41l6XxgYEB9Pf3l2w3MjJi1LjIavjdu3djeHi4ZErZxWbBkyV5uZ25Z755CJwcC5/L5UoyZfl8HoqiIBKJYGJiwnJfRKNR5HI5x6YaIloGgtY1XdcFAJHJZIQQQmQyGQFAJJNJ2+3lejNVVYWiKI7rzccvlwZVVYWqqmXTC0Doui40TTPeUwghFEURyWRSaJpmOUYymbRsV6mV/mqYr1vxZyIVXx+5na7rQojS66yqqshkMsax5b+RSMT4v/mlaVpJesz7mddrmma8l/l9zfeOpmnGOamquqLXtKenR/zgBz9YsfcjWq0Y5MkiEomUDYrmYJRMJo2gIV/mgCCDRLnjyUAlA0i595cBW+5jzojIv4sDXXFmoFKNCvLF16OYoihGoC8O+k6ZqeKMj10mTlVV2/eMRCLGZ+ykXCZO7r9Y5q3eGOSJFjDIrzOKopSU4Cp9JZNJSwmvOHjK0nQlAVIGM3NJVDK/h3mdzEREIhGhaZrlGDL4yUBVnO5qNSrIlwvwkqIotkFXBlS7V7maFCGcg7w5bU6B2u745kxeIzDIEy24rur6fVrTnDqrVUq28+q6XtLRLpfLlR0uJdtrARgTvNiR86Kbtwc+nCXP3IbtdAyn5auVPOdKOH2GiURiWdq7Zdpkm79TXw3z5yU78RFRY3kEv4lEJTweD4PUGrZr1y48/PDD2LlzZ6OTQtRQ7F1PRETkUgzyRERELsUgT0RE5FIM8kRERC7FIE9ERORSDPJEREQuxSBPRK7DX6EjWsAgT0RE5FIM8kRERC7FIE9ERORSDPJEREQuxSBPRETkUgzyRERELsUgT0SuwyF0RAsY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjIpRjkiYiIXIpBnohch73riRYwyBMREbkUgzwREZFLMcgTERG5FIM8ERGRSzHIExERuRSDPBERkUsxyBMREbkUgzwRuQ7HyRMtYJAnIiJyKQZ5IiIil2KQJyIicikGeSIiIpdikCci19myZUujk0C0KniEEKLRiSBabTweD/jVIKK1jiV5Ihv33HNPo5NARLRkLMkTERG5FEvyRERELsUgT0RE5FIM8kS0rqRSKXg8HqRSqYq39/v9Jcuj0Sii0Sjy+Tw8Hk/F7+/3+5HNZivaNhQKIRqN2q7zeDyLHicejyMUCpUsj0ajFaXZ7/cjHo/bLq/0+lFjMciTa9g99FKplONDshLZbLaih6EMHMXBo9L93Uhej1pelQbBauXzeXR3dwMA+vv763LMK1euQFVVeDwe5PP5Rd9f1/WKjpvP55FOp5FIJErWpVIpqKqKQCBQ9hinTp1CR0eHZVk2m8Xg4CAURSn73ZBp3b17d0XppdWJQZ5WLXPgNL+qsXXrVkxMTBilmXg8vmiAsSu5SLLUVvwaGRmBEAJCCGzdurXq/eXLrtS1VsnrUfxSFAWRSMRxvRBi0eBVq87OTqiqagyPdPqso9GobendTiAQwPj4OCKRCI4cOWLsb/f5KooCAAgGg7br5XuGQiFjW7vMT39/P9LpdNl7V2YSzEE6n88jGAwimUwil8thcHDQsUQ+OTkJVVUxOTlZ8j66rqO7u7tkOUv3qw+DPK1qiqIYD/5MJlP1/j6fD7lcDsBCpqG3t7dscFFVtaLj6rpu2W98fNzyntXuL4SApmlVnx9VTmag5Gc1MTGBWCxmW2uQSCSg63pVQSuRSBil7kQi4Xh/qarqmMmR9yoAaJpmWSfJ5gO7zJPZ6OgoVFU17sd8Pg9FUZBMJhEOhwHACNZ212BkZAQ9PT0Ih8O275VMJkuWy+PS6sEgT2uOU2kc+LCEFI/HLe2R4+PjVT2AZDV7MBgEAMt7VJPG4v2vXLlS8TGofkKhENLptCWI+nw+ZDIZBINB2yAXiUQsVfrFpfPBwUEMDg5W3NQgq8bHx8exefPmmpuR+vv7K2pqGB4eNgJ8PB6HoijIZDKW74HP54Ou6wgGg5b0yFqAcrVStDYwyNOa41QaB4BMJgMhBHp7e0v28/v9ZatIzQKBgKX2oLg0VWkai/ffuHFjLafsOhMTEytWvevxeJDL5YzPz1wVHwgEjEBfnJb9+/dD13VcunQJQGnpXJbKK2lqkJkC2T7e29uLzZs3V9wkIB0/fhx+vx89PT2262+++WYAC6V92fYfjUYxPDzsmDafzwchhKVZSzY7SMUZa7vq+uXqR0FLwyBPq5qu68ZDRJaKgYWSWS0BoriKMZlM1jO5VCG/349MJoPu7u5l64cg+3RommYpwReTGbru7m5L0JXBTwbOYrlcDul0umwaZMZSZj5jsZhRYu7t7TUyO5UG+wcffBDj4+Po7Owse936+/uNpqdEImGcf3GvePPIgVwuh/HxceTzeQwODpYc0y5DU23ml1Yegzytak5t8j09PXXrHV2JaDRqeagqilK2A10xWRJidf0Cn89nBNdcLlfSM92p02Wlr2984xvo7u52rNWxI4TA3r17LWkx98Y3kz3PVVW1zWzKdExMTFhK0LLELPcxl6LlPvK9Y7GY4z0lA7ddx8FUKoXOzs6K+obYOXLkCCKRSMlyu45+tXSGpZXFIE9rkmxXrLaKsLiK0e4BLmWzWaP2YPv27ZbOdeZSjKwWNXemAxYetubaB1bX28vlctA0DYqiGEHLrrNXNa9vfetbNZUwZTOLz+ezdFQrJju19fT0YGRkpGS9TIddoM3lciX9Q2SwN+/j1PFOGhkZQSwWs2SONm3ahDNnzmD//v1Vn7uUz+dt92dJfm1ikKc1a+/evRgYGKhqn0qr61OpFAYGBozaA/lQrrQkHo1GcebMGWP/SkuT61Vvby90XXfs7b7SZIB3GtM+PDxs9DxPp9O2aY7H4yXV8NFotG79EAKBACKRCEZHRwHASGsikShbit+0aVPZ45ozs2Ysya9NDPK0ZvX29i7aJlqrcDhs+7C7fPlyyVAlO+bhVFQZWZpdrjHy1ejs7EQmk7ENlubaBmChxG2X2RweHsbevXuNuRFkCbm7u7ukmr2aWfDMEolExRlI2cRQa20SS/JrE4M8rWkr/YA5c+YMOjs7V/Q9aeXlcjnbzEY2m0UsFsPExISxrLe3F7lczjIELZvNQtf1kgAsMzLDw8PI5/NG34ORkREEg8GyEzFVolwAn5ychKIoNbfVsyS/Nl3X6AQQlSN719dLd3d3STt8JSVzaXBwsKZJeWjtk300kslkSaCUnQeBhdJ1T0+Pbec18/bAh53cAoEAdF2HoijYtm0bgIWOd7FYrKK0mTvyOenv78fevXsrOp4dVVUdq/IZ6FcvBnla1RRFMR6I5o5wtTLP9lWtaDQKRVGqrk62y0Q4ZSwqnXGPVpacMKbc/SOEgMfjwebNm0sm3VFV1fEzl23pspQvaZpmqQkwB9J4PF6SASg3HDQajdrWLCzmypUrlnSXC+bBYLBsRoAag0GeVq1wOGx5oMohV1Lxg26xH+woN1a6mOx4BcAokQ0ODlb84yLmtNmV6HRdLyl1xeNxnDp1quI00soovg/LcWo+qjbw2W1vPnZvb29VbfHV3LtmGzduZJv7GucR/ATJpWTQLPeAtasd4FeCiNyCQZ6IiMil2LueiIjIpRjkiYiIXIpBnoiIyKUY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjIpRjkiYiIXIpBnoiIyKUY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjIpRjkiYiIXIpBnoiIyKUY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjIpRjkiYiIXIpBnoiIyKUY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjIpRjkiYiIXIpBnoiIyKUY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjI5fPu4wAAAEJJREFUpRjkiYiIXIpBnoiIyKUY5ImIiFyKQZ6IiMilGOSJiIhcikGeiIjIpRjkiYiIXIpBnoiIyKUY5ImIiFzq/wNDHMf27z2dvgAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 3,
     "status": "ok",
     "timestamp": 1649954662902,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "8sywqMFs9fWd",
    "outputId": "d5c626fd-70c9-44f7-a4c3-3b112e4654c4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MRP中每个状态价值分别为\n",
      " [[-2.01950168]\n",
      " [-2.21451846]\n",
      " [ 1.16142785]\n",
      " [10.53809283]\n",
      " [ 3.58728554]\n",
      " [ 0.        ]]\n"
     ]
    }
   ],
   "source": [
    "def compute(P, rewards, gamma, states_num):\n",
    "    ''' 利用贝尔曼方程的矩阵形式计算解析解,states_num是MRP的状态数 '''\n",
    "    rewards = np.array(rewards).reshape((-1, 1))  #将rewards写成列向量形式\n",
    "    value = np.dot(np.linalg.inv(np.eye(states_num, states_num) - gamma * P),rewards)\n",
    "    return value\n",
    "\n",
    "\n",
    "V = compute(P, rewards, gamma, 6)\n",
    "print(\"MRP中每个状态价值分别为\\n\", V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 340,
     "status": "ok",
     "timestamp": 1649954667427,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "5ILxWaLR9fWd"
   },
   "outputs": [],
   "source": [
    "S = [\"s1\", \"s2\", \"s3\", \"s4\", \"s5\"]  # 状态集合\n",
    "A = [\"保持s1\", \"前往s1\", \"前往s2\", \"前往s3\", \"前往s4\", \"前往s5\", \"概率前往\"]  # 动作集合\n",
    "# 状态转移函数\n",
    "P = {\n",
    "    \"s1-保持s1-s1\": 1.0,\n",
    "    \"s1-前往s2-s2\": 1.0,\n",
    "    \"s2-前往s1-s1\": 1.0,\n",
    "    \"s2-前往s3-s3\": 1.0,\n",
    "    \"s3-前往s4-s4\": 1.0,\n",
    "    \"s3-前往s5-s5\": 1.0,\n",
    "    \"s4-前往s5-s5\": 1.0,\n",
    "    \"s4-概率前往-s2\": 0.2,\n",
    "    \"s4-概率前往-s3\": 0.4,\n",
    "    \"s4-概率前往-s4\": 0.4,\n",
    "}\n",
    "# 奖励函数\n",
    "R = {\n",
    "    \"s1-保持s1\": -1,\n",
    "    \"s1-前往s2\": 0,\n",
    "    \"s2-前往s1\": -1,\n",
    "    \"s2-前往s3\": -2,\n",
    "    \"s3-前往s4\": -2,\n",
    "    \"s3-前往s5\": 0,\n",
    "    \"s4-前往s5\": 10,\n",
    "    \"s4-概率前往\": 1,\n",
    "}\n",
    "gamma = 0.5  # 折扣因子\n",
    "MDP = (S, A, P, R, gamma)\n",
    "\n",
    "# 策略1,随机策略\n",
    "Pi_1 = {\n",
    "    \"s1-保持s1\": 0.5,\n",
    "    \"s1-前往s2\": 0.5,\n",
    "    \"s2-前往s1\": 0.5,\n",
    "    \"s2-前往s3\": 0.5,\n",
    "    \"s3-前往s4\": 0.5,\n",
    "    \"s3-前往s5\": 0.5,\n",
    "    \"s4-前往s5\": 0.5,\n",
    "    \"s4-概率前往\": 0.5,\n",
    "}\n",
    "# 策略2\n",
    "Pi_2 = {\n",
    "    \"s1-保持s1\": 0.6,\n",
    "    \"s1-前往s2\": 0.4,\n",
    "    \"s2-前往s1\": 0.3,\n",
    "    \"s2-前往s3\": 0.7,\n",
    "    \"s3-前往s4\": 0.5,\n",
    "    \"s3-前往s5\": 0.5,\n",
    "    \"s4-前往s5\": 0.1,\n",
    "    \"s4-概率前往\": 0.9,\n",
    "}\n",
    "\n",
    "\n",
    "# 把输入的两个字符串通过“-”连接,便于使用上述定义的P、R变量\n",
    "def join(str1, str2):\n",
    "    return str1 + '-' + str2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 3,
     "status": "ok",
     "timestamp": 1649954670178,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "juDFPGkP9fWe",
    "outputId": "20903c97-0f0e-4fb2-93e1-5b6b650fda1a"
   },
   "outputs": [],
   "source": [
    "gamma = 0.5\n",
    "# 转化后的MRP的状态转移矩阵\n",
    "P_from_mdp_to_mrp = [\n",
    "    [0.5, 0.5, 0.0, 0.0, 0.0],\n",
    "    [0.5, 0.0, 0.5, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0, 0.5, 0.5],\n",
    "    [0.0, 0.1, 0.2, 0.2, 0.5],\n",
    "    [0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "]\n",
    "P_from_mdp_to_mrp = np.array(P_from_mdp_to_mrp)\n",
    "R_from_mdp_to_mrp = [-0.5, -1.5, -1.0, 5.5, 0]\n",
    "\n",
    "V = compute(P_from_mdp_to_mrp, R_from_mdp_to_mrp, gamma, 5)\n",
    "print(\"MDP中每个状态价值分别为\\n\", V)\n",
    "\n",
    "# MDP中每个状态价值分别为\n",
    "#  [[-1.22555411]\n",
    "#  [-1.67666232]\n",
    "#  [ 0.51890482]\n",
    "#  [ 6.0756193 ]\n",
    "#  [ 0.        ]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 317,
     "status": "ok",
     "timestamp": 1649954714601,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "3gKVFNen9scC",
    "outputId": "3d5b5f1b-d5d2-4a26-fde2-76843910507b"
   },
   "outputs": [],
   "source": [
    "def sample(MDP, Pi, timestep_max, number):\n",
    "    ''' 采样函数,策略Pi,限制最长时间步timestep_max,总共采样序列数number '''\n",
    "    S, A, P, R, gamma = MDP\n",
    "    episodes = []\n",
    "    for _ in range(number):\n",
    "        episode = []\n",
    "        timestep = 0\n",
    "        s = S[np.random.randint(4)]  # 随机选择一个除s5以外的状态s作为起点\n",
    "        # 当前状态为终止状态或者时间步太长时,一次采样结束\n",
    "        while s != \"s5\" and timestep <= timestep_max:\n",
    "            timestep += 1\n",
    "            rand, temp = np.random.rand(), 0\n",
    "            # 在状态s下根据策略选择动作\n",
    "            for a_opt in A:\n",
    "                temp += Pi.get(join(s, a_opt), 0)\n",
    "                if temp > rand:\n",
    "                    a = a_opt\n",
    "                    r = R.get(join(s, a), 0)\n",
    "                    break\n",
    "            rand, temp = np.random.rand(), 0\n",
    "            # 根据状态转移概率得到下一个状态s_next\n",
    "            for s_opt in S:\n",
    "                temp += P.get(join(join(s, a), s_opt), 0)\n",
    "                if temp > rand:\n",
    "                    s_next = s_opt\n",
    "                    break\n",
    "            episode.append((s, a, r, s_next))  # 把（s,a,r,s_next）元组放入序列中\n",
    "            s = s_next  # s_next变成当前状态,开始接下来的循环\n",
    "        episodes.append(episode)\n",
    "    return episodes\n",
    "\n",
    "\n",
    "# 采样5次,每个序列最长不超过1000步\n",
    "episodes = sample(MDP, Pi_1, 20, 5)\n",
    "print('第一条序列\\n', episodes[0])\n",
    "print('第二条序列\\n', episodes[1])\n",
    "print('第五条序列\\n', episodes[4])\n",
    "\n",
    "# 第一条序列\n",
    "#  [('s1', '前往s2', 0, 's2'), ('s2', '前往s3', -2, 's3'), ('s3', '前往s5', 0, 's5')]\n",
    "# 第二条序列\n",
    "#  [('s4', '概率前往', 1, 's4'), ('s4', '前往s5', 10, 's5')]\n",
    "# 第五条序列\n",
    "#  [('s2', '前往s3', -2, 's3'), ('s3', '前往s4', -2, 's4'), ('s4', '前往s5', 10, 's5')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 292,
     "status": "ok",
     "timestamp": 1649954717890,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "uZR44aSO9fWf",
    "outputId": "7354c75f-1bc7-44b2-accc-85019278720d"
   },
   "outputs": [],
   "source": [
    "# 对所有采样序列计算所有状态的价值\n",
    "def MC(episodes, V, N, gamma):\n",
    "    for episode in episodes:\n",
    "        G = 0\n",
    "        for i in range(len(episode) - 1, -1, -1):  #一个序列从后往前计算\n",
    "            (s, a, r, s_next) = episode[i]\n",
    "            G = r + gamma * G\n",
    "            N[s] = N[s] + 1\n",
    "            V[s] = V[s] + (G - V[s]) / N[s]\n",
    "\n",
    "\n",
    "timestep_max = 20\n",
    "# 采样1000次,可以自行修改\n",
    "episodes = sample(MDP, Pi_1, timestep_max, 1000)\n",
    "gamma = 0.5\n",
    "V = {\"s1\": 0, \"s2\": 0, \"s3\": 0, \"s4\": 0, \"s5\": 0}\n",
    "N = {\"s1\": 0, \"s2\": 0, \"s3\": 0, \"s4\": 0, \"s5\": 0}\n",
    "MC(episodes, V, N, gamma)\n",
    "print(\"使用蒙特卡洛方法计算MDP的状态价值为\\n\", V)\n",
    "\n",
    "# 使用蒙特卡洛方法计算MDP的状态价值为\n",
    "#  {'s1': -1.228923788722258, 's2': -1.6955696284402704, 's3': 0.4823809701532294,\n",
    "# 's4': 5.967514743019431, 's5': 0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 303,
     "status": "ok",
     "timestamp": 1649954723228,
     "user": {
      "displayName": "Sam Lu",
      "userId": "15789059763790170725"
     },
     "user_tz": -480
    },
    "id": "COkP4ZDh9fWg",
    "outputId": "943a07c9-f8db-4646-841b-d2960785eb17"
   },
   "outputs": [],
   "source": [
    "def occupancy(episodes, s, a, timestep_max, gamma):\n",
    "    ''' 计算状态动作对（s,a）出现的频率,以此来估算策略的占用度量 '''\n",
    "    rho = 0\n",
    "    total_times = np.zeros(timestep_max)  # 记录每个时间步t各被经历过几次\n",
    "    occur_times = np.zeros(timestep_max)  # 记录(s_t,a_t)=(s,a)的次数\n",
    "    for episode in episodes:\n",
    "        for i in range(len(episode)):\n",
    "            (s_opt, a_opt, r, s_next) = episode[i]\n",
    "            total_times[i] += 1\n",
    "            if s == s_opt and a == a_opt:\n",
    "                occur_times[i] += 1\n",
    "    for i in reversed(range(timestep_max)):\n",
    "        if total_times[i]:\n",
    "            rho += gamma**i * occur_times[i] / total_times[i]\n",
    "    return (1 - gamma) * rho\n",
    "\n",
    "\n",
    "gamma = 0.5\n",
    "timestep_max = 1000\n",
    "\n",
    "episodes_1 = sample(MDP, Pi_1, timestep_max, 1000)\n",
    "episodes_2 = sample(MDP, Pi_2, timestep_max, 1000)\n",
    "rho_1 = occupancy(episodes_1, \"s4\", \"概率前往\", timestep_max, gamma)\n",
    "rho_2 = occupancy(episodes_2, \"s4\", \"概率前往\", timestep_max, gamma)\n",
    "print(rho_1, rho_2)\n",
    "\n",
    "# 0.112567796310472 0.23199480615618912"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "第3章-马尔可夫决策过程.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
